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Weakly Interactive Massive Particles (WIMPs) are the most widely studied candidate particles 
forming the cold dark matter (CDM) whose existence can be inferred from a wealth of astrophys- 
ical and cosmological observations. In the framework of the minimal cosmological model detailed 
measurements on the cosmic microwave background by the PLANCK collaboration fix the scaled 
CDM relic density to Q,ch? = 0.1193 ± 0.0014, with an error of less than 1.5%. In order to fully 
exploit this observational precision, theoretical calculations should have a comparable or smaller er¬ 
ror. In this paper we use recent lattice QCD calculations to improve the description of the thermal 
plasma. This affects the predicted relic density of “thermal WIMPs”, which once were in chemical 
equilibrium with Standard Model particles. For WIMP masses between 3 and 15 GeV, where QCD 
effects are most important, our predictions differ from earlier results by up to 9% (12%) for pure 
S'—wave (P—wave) annihilation. We use these results to compute the thermally averaged WIMP 
annihilation cross section that reproduces the correct CDM relic density, for WIMP masses between 
0.1 GeV and 10 TeV. 


1. INTRODUCTION 

Assuming Newtonian gravity, and its extension into General Relativity, describes gravitational forces correctly, 
astronomical and cosmological observations show that most of the matter in our Universe is a non-luminous, neutral 
form of matter, called dark matter (DM). Quantitatively, it accounts for ^ 85% of all matter [IHS]. There are many 
possible DM candidates [1]. Among those, weakly interactive massive particles (WIMPs) have been studied in most 
detail. There are several reasons for the popularity of WIMPs. If they once were in full chemical equilibrium with the 
particles of the Standard Model (SM), which is true if the largest temperature after the most recent period of entropy 
production exceeded about 5% of the WIMP mass, the WIMP relic density can be calculated unambiguously for a given 
cosmological model, independent of initial conditions, using only particle physics quantities (masses and couplings, or 
cross sections) as input. This calculation yields approximately the observed relic density for roughly weak-strength 
WIMP annihilation cross sections. This not only hints at a deep connection between DM and extensions of the SM 
addressing some of its shortcomings [5H9], it also allows various ways in which the existence of WIMP DM can be 
probed. Unfortunately these probes have so far not yielded an unambiguous signal. Successful WIMP candidates 
must therefore not only satisfy the relic abundance constraint, but also constraints coming from direct unHH] and 
indirect [HE] detection experiments. 

In standard cosmology it is assumed that the comoving entropy density remained constant since the epoch when 
WIMPs were in full thermal equilibrium with SM particles. The quantity of interest is then the ratio of the WIMP 
number density and the entropy density s. In order to compute the temperature dependence of n^, one thus needs 
to know the precise temperature dependence of s. In addition, the expansion rate, described by the Hubble parameter 
H, is proportional to the square root of the total energy density p; hence the temperature dependence of p also has 
to be known. In early analyses OE] s{T) and p{T) were calculated ignoring all interactions between SM particles, 
i.e. treating the thermal plasma as a relativistic free gas. This is quite a good approximation for most temperatures. 
However, it was realized early on [181 II9j that this approach fails for temperatures near the deconfinement transition, 
i.e. the transition from hadronic (pions, kaons, ...) to partonic (quarks, gluons) degrees of freedom. More recently 
it was realized that also for some range of temperatures above this transition the interactions between quarks and 
gluons should not be ignored HOI mi- 

In this paper we carefully model the effect of strong interactions on the temperature dependence of the energy 
and entropy densities. We exploit recent calculations performed in the framework of the Hadron Resonance Gas 
(HRG) model well below the deconfinement transition [52], smoothly matching to results from lattice QCD (LQCD) 
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at T = 100 MeV [531. that this can change the predicted relic density by more than 5% compared to earlier 

treatments [l9l - l2T] of strong interaction effects on the thermal plasma. 

We then use our improved treatment of strong interaction effects to determine the value of (cru), the thermally 
averaged product of WIMP annihilation cross section times relative velocity of the annihilating WIMPs, that is 
required to reproduce the correct relic density, under the assumption that (cru) is independent of temperature. As 
already noted in ref. [21], depending on the WIMP mass, this product can deviate by up to a factor of ^ 1.5 in either 
direction from the “canonical” value of 3 • 10“^® cm®s“^. This is important since experiments are now beginning to 
be sensitive to the canonical cross section. For example, analyses of FERMI-LAT gamma ray observations of nearby 
dwarf galaxies and of the diffuse gamma ray emission in the galaxy have, for some combinations of WIMP masses 
and WIMP annihilation final states, constrained the WIMP annihilation cross section to be lower than the canonical 
one [SHEj. 

The remainder of this paper is organized as follows. In Sec. [^ we review the calculation of the relic abundance 
of thermal WIMPs. We outline earlier estimates of s(T) and p{T), and explain in detail how we compute these 
quantities. In Sec. we present our result for (av), and compare with previous results. In Sec. we compare our 
{av) with bounds on this quantity from indirect WIMP detection experiments and CMB anisotropies. Finally, we 
summarize the main results of our work in Sec. [5| 


2. CALCULATION OF THE RELIC ABUNDANCE 
2.1. Basic Pramework 


The starting point of our calculation is the Boltzmann equation, which describes the time evolution of the number 
density of a stable particle species x that can annihilate into a pair of SM particles [18] : 


dny , , 

— +3ilnx = -{av) 


(2 2 
("x - ^X.eq 
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( 1 ) 


Here t is the cosmological time, H is the Hubble parameter, (av) is the thermal average of the product of WIMP 
annihilation cross section and Mpller velocity^, and is the density of x particles in thermal equilibrium. We 

have assumed that x particles are self-conjugate (Majorana) particles. With minor modifications our results can also 
be applied to particles that are not self-conjugate, such as Dirac fermions or complex scalar particle. In this case one 
will in general have to track the densities of particles and antiparticles separately. 

Since the total entropy is assumed to be conserved, S = sa^ =const. where a is the scale factor in the Friedman- 
Robertson-Walker metric and s is the entropy density, it is convenient to define = n^fs. Using the definition 
H = {l/a)da/dt it is easy to see that this definition absorbs the dilution term 3Hn^ on the right-hand side of the 
Boltzmann eq.Q. 

Moreover, since s as well as n^^eq depend explicitly on temperature rather than on time, one replaces time t by 
X = m^/T. To that end, the entropy density is written as 


27r2 

<T) = -^h{T)T^ . 


( 2 ) 


Here h{T) counts the effective number of relativistic degrees of freedom (d.o.f.) that contribute to s, i.e. if all relevant 
particles have a common temperature T and masses ^ T, h{T) —>■ gi, where gi is the number of effective internal 
degrees of freedom for particle species i. For example, a massless photon contributes g.y = 2, a massless Dirac fermion 
(with both helicities being in thermal equilibrium) contributes ge = 4 • 7/8 = 3.5, and so on. Using ^ [/i(T)T®a®] = 0 
one then finds 


dT _ TH 

dt 1 -I- I ^ inh(T) 

^ 3 d In T 


(3) 


We thus also need to know the temperature dependence of the Hubble parameter. To that end we use the Friedmann 
equation = SttGnp/^, where p and are the energy density and Newton’s gravitational constant, respectively. 


^ In the non-relativistic limit the M0ller velocity becomes the relative velocity between the annihilating WIMPs. 
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During the radiation dominated epoch, in which the decoupling of WIMPs falls, the energy density can be written as 

P{T) = ^^9iT)T^ ■ (4) 

Here g{T) counts the effective number of relativistic d.o.f. that contribute to p. In the limit where all particle masses 
can be ignored, g{T) = h{T), i.e. the numbers of d.o.f. defined via the entropy density and via the energy density are 
the same, but in general this is not the case. 

Putting everything together, we have 


where we have introduced 


and 
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( 6 ) 

(7) 


being the WIMP mass. Numerically, A = 2.76 x 10® for a WIMP mass of 1 GeV and (av) = 10“^® cm^s“^. The 
quantities g and h appearing in eq.([^ are the effective number of relativistic d.o.f. defined via the energy density, 
eq.Q, and entropy density, eq.([^, respectively. Note that we need to know both g and h: both appear in the 
definition of gV^^ and h also appears in the explicit expression for the scaled equilibrium density WIMPs 

decouple when they are essentially non-relativistic, so we can write 


Y^^eq{x) — 
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where = 2 for a neutral Majorana fermion. 


2.2. The Functions g{T) and h{T) 


Most papers on the calculation of the WIMP relic density focus on the annihilation cross section. Our concern 
is instead an accurate calculation of the effective numbers of degrees of freedom encoded in the functions g(T) and 
h{T). To that end, we have to compute the energy density p{T) and the pressure p{T); the entropy density is 
then given by s{T) = [p{T) +p{T)]/T. If the single particle distribution functions fi{k,T) for a particle species i 
is known, the contribution of this species to the energy density and pressure are given by HH] PiiT) = gi/(27r)3 . 
f d^kEi{k)fi{k,T), Pi{T) = gi/(27r)^ • / d^k{k)'^f{k,T)/{3Ei{k)); here gi is again the number of degrees of freedom 
of species i, and k is the three-momentum. 

The simplest case is that of a free (non-interacting) particle with mass rrii. The distribution function is then the 

Bose-Einstein or Fermi-Dirac distribution function, which depends only on the ratio of the energy Ei = mf + 
and the temperature. Introducing the dimensionless quantities yt = Eijrrii and Xi =mi/T, one has [T^fTO] : 
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( 10 ) 


In the denominators of the integrals +1 applies to fermions and —1 to bosons. 

In the presence of strong interactions, Eqs.(|^ and (10) no longer provide good approximations. Before describing 
our own calculation of these functions, we briefly review the state of the art, as encoded in the widely used program 
packages for the calculation of the WIMP relic density DarkSUSY [26], micrOMEGAs [27] and Superlso [28] . 

It was realized quite early on that care has to be taken to describe the thermodynamics of the early universe around 
the deconfinement transition. In [29] the interactions between hadrons and between partons were approximated by 
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simple non-relativistic potentials. Ref. |30j instead used free particles, and defined the transition temperature from 
the hadronic to the partonic phase as that temperature where the two calculations give the same entropy density; 
note that the hadronic phase gives a quickly rising h(T), since the number of contributing hadrons quickly increases 
with temperature. This “hadron resonance gas model” was used in all subsequent calculations at sufficiently low 
temperatures, including our own. 

One problem of the simple definition of the transition temperature used in ref. |30j is that it leads to a discontinuity 
in g{T). Ref.p5] therefore used smooth functions interpolating between the hadronic and partonic phases. While 
these functions ensure that not only g{T) and h{T), but also their derivatives are smooth, they were not based on 
dynamical considerations. The authors advocated estimating the uncertainty by using two quite different values, 150 
and 400 MeV, for the transition temperature. The same functions were used in Ref. [19] : the functions for a transition 
temperature of 150 MeV are still used by default in the computer packages mentioned above. 

The hrst attempt to include the results of lattice QCD calculations was due to Hindmarsh and Philipsen [5D] . At 
that time the most accurate lattice QCD calculations did not include dynamical quarks. There was some evidence 
that the ratio of the true pressure to the corresponding value for non-interacting particles shows little dependence on 
the number of quark flavors |3I]. Ref. [20] therefore scaled the contribution of all strongly interacting partons by the 
same correction function, determined from pure glue lattice calculations |3T]; at T = 1.2 GeV, these were matched to 
perturbative calculations |32j . 

The treatment by Laine and Schroeder |21j is rather similar. However, their results are based on a different set 
of pure glue lattice QCD calculations |53|. Moreover, they match to perturbative calculations at the much lower 
temperature of 350 MeV. Finally, they include the quark mass dependence up to next-to-leading order, 0{g'^), in the 
perturbative expansion. In particular, they point out that charm quarks make non-negligible contributions already 
at temperatures of a few hundred MeV. 

We now describe own treatment. At temperatures well above the electron mass, we treat all SM particles without 
strong interactions as free particles, i.e. we use Eqs.([^ and (10) for these particles. This includes the leptons, the 
electroweak gauge bosons as well as the single physical Higgs boson of the SM. Note that for the physical Higgs 
mass, rriff ~ 125 GeV, in the Standard Model electroweak symmetry breaking leads to a smooth cross-over, not 
a phase transition |34H36j : hence the comoving entropy density remains constant, as assumed in the derivation of 
eq. For simplicity we compute g{T) and h{T) using free, massive and Z bosons (with three d.o.f. each) and a 
single physical Higgs boson even for temperatures above the electroweak cross-over, where a more accurate treatment 
would use massless gauge bosons (with two d.o.f. each) and a massive complex Higgs doublet (with four d.o.f.). The 
difference between these treatments is very small, and will only affect the relic density of WIMPs with masses above 
2 TeV. 

The main focus of our work is on the effect of QCD interactions. These are most important around the deconfinement 
transition, which was also a smooth crossover m with conserved total entropy. We use the results of a recent lattice 
calculation with Nf =2 + 1 active flavors (meaning equal masses are used for u and d quarks, but the larger mass of 
the strange quark has been taken into account) [23]. This calculation covers temperatures between 100 and 400 MeV. 
Ref. [23] provides a parameterization of the pressure due to u, d, s quarks and gluons in this temperature range: 


^ ^ [1 
Ti 2 ^ 


tanh (ct(f - to))] ’ 


Pid + an/t + bn/P + dn/t^ 
1 + O'd/t + bd/P + dd(P 


( 11 ) 


Here i = T/Tc, Tc = 154 MeV being the QCD transition temperature. In this parameterization, pid = 197 r ^/36 is 
the ideal gas value of p/T'^ fo r QC D with three massless quarks. The values of the numerical coefficients appearing 
in eq.(ll) are listed in Table 2.2 Recall that tanh(z) approaches unity for large argument z. Therefore eq.(ll) 
automatically approaches the ideal gas value for T :+ T^, i.e. t :+ 1. Moreover, ref.|23| shows that eq. 0 matches 
quite well to available perturbative calculations at higher temperature. We therefore use this parameterization to 
describe the contribution from u, d, s quarks and gluons for all temperatures above 100 MeV. 


Ct 


bn 

dn 

3.8706 

-8.7704 

3.9200 

0.3419 

to 

ad 

bd 

dd 

0.9761 

-1.2600 

0.8425 

-0.0475 


TABLE I: Parameters used in eq. (Ill to describe the pressure of (2+l)-fiavor QCD. 


Once the pressure p is known, the energy density p can be computed from the relation between the trace of the 
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energy-momentum tensor, also called the trace anomaly, and the pressure 


I{T) p-3p 




'P'i 


A 

dT 




( 12 ) 


Since we have the analytical expression (11) for pjT^, we can easily obtain its derivative, and thus I(T); the first 
equation in ( |I^ then allows to compute p(T) from I{T) and p{T). Finally, once p{T) and p{T) are known, s(T) = 
[p{T) + p{T)UT can also be computed. 

As noted above, the effect of the charm quark is not negligible at temperatures near [H]. We included its 
contribution to the functions g and h using lattice QCD results from Table 6 of [SS], using the physical ratio of 
charm and strange quark masses, rric/ras = 11.85 [39]. This gives us pdT'^, from which the charm contributions to 
p and s can be obtained as outlined in the previous paraCTaph. Thi s description is valid for T < 1 GeV. For larger 
temperatures we smoothly match to the ideal gas results and ( 10 ), using a fit function like ( 11 ) with pid = jQQ 


and different values for coefficients and powers to interpolate between the two regimes. This ensures that not only 
the functions g and h but also their first derivatives are smooth everywhere. 

Bottom and top quarks contribute significantly only at high temperatures, where even QCD interactions have 
become relatively small. We therefore treat these quarks as free particles, with on-shell masses given by the Particle 
Data Collaboration [4(1] . 

At temperatures lower than Tc, the thermodynamic behavior of QCD can be described by the hadron resonance gas 
model, in which all the hadrons and hadron resonances are considered to contribute to the thermodynamics quantities 
as non-interacting particles. As noted above, this has been used already in the early treatments [211 ISO]- Ref. [23] 


shows that for temperatures between 100 MeV and Tc it matches well to the QCD results parameterized in eq.(ll|. 
A convenient parameterization of the trace anomaly in this model can be found in ref. [ 22 ] : 


I{T) _P ^ ^ 






(13) 


with oi = 4.654 GeV“^, 02 = —879 GeV“^, 03 = 8081 GeV”"^, 04 = —7039000 GeV“^°. This parameterization is valid 
for 70 MeV < T < Tc. We use it to describe the contribution from strongly interacting particles for all temperatures 
T < 100 MeV, using cubic splines to interpolate smoothly to QCD results at T > 100 MeV. The contribution from 
charmed particles is negligible at these low temperatures.^ 

By inverting eq. (12), the pressure can be calculated: 


P{T) _ Po 


n 


'To 


pi5 


(14) 


using the numerical result p(To)/Tq = 0.1661 at Tq = 70 MeV [22]. The integral in eq.((l4]) can easily be evaluated 
analytically if I{T) is given by eq.(13). Eq.(14| thus again provides us with an analytical parameterization of the 


pressure, from which we can compute the energy and entropy densities as described above. 

At temperatures below 1 MeV, the effect of neutrino decoupling should be included. The rate of reactions changing 
the and Vt number densities actually becomes smaller than the Hubble parameter, indicating decoupling, at a 
temperature of several MeV. However, at first the expansion of the universe affects photons and neutrinos in the same 
way even after neutrino decoupling, i.e. the photon and neutrino temperatures remain the same. This changes only 
once e“'"e“ pairs begin to annihilate, at T ~ me- Since neutrinos are already (almost) decoupled by this time, the 
entropy that was stored in electrons and positrons gets transferred (almost) entirely to photons, not to neutrinos. 
In the limit where neutrino decoupling was complete when electron decoupling began, this argument shows that for 
T <C rUe the ratio of relic photon and neutrino temperatures is Ty/T^ = (11/4)^/^. Actually (electron) neutrinos were 
not completely decoupled at T ~ rUe- This can be described by writing 


h = 2 


7 [A 

1+o^eff ( YY 
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(15) 


g = 2 
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(16) 


^ Strictly speaking the hadronic contribution to g and h should become exponentially small at T <C m-Tr = 140 MeV. At very low 
temperatures eq. therefore is not accurate. However, this is not important for us, since for T <§C 100 MeV the hadronic contribution 
is in any case very small; this small contribution need not be described very accurately. 
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with iVeff ~ 3.046. Note that these expressions include the contribution from the photon, with gy = 2. Eqs.(15) and 
(16) are applicable for T <C We, in practice for T < 50 keV. As mentioned above, for T > 1 MeV we have = T^. 
Fot 50 keV <T < 1 MeV we use numerical results from Fig. 1 of m to determine the evolution of T^, with respect 
to T^. This can then be plugged into eqs.(15) and (16) instead of T^/T^ = (4/11)^/^ to compute the photon and 


neutrino contribution to g{T) and h{T). Note that the temperature T is defined to be that of the photons, T = T. 


7- 




FIG. 1: The functions h{T) (top frame) and gV^{T) (bottom) defined in eqs. ([^ and The original calculation by Gondolo 
and Gelmini m, based on results from ref. [^, are shown by the green dashed curves. The red dot-dashed and blue dotted 
curves show results from refs. and m, which are based on pure glue lattice QCD calculations. The black solid curves 
depict our results, which are based on lattice calculations with Nf = 2 + 1 dynamical quark flavors. 
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3. RESULTS AND COMPARISON WITH PREVIOUS STUDIES 


We are now ready to present some numerical results. Figure shows the functions h{T) and gl^^(T) that parame¬ 
terize thermodynamic effects in the Boltzmann equation ([^. The black solid curves show our results, while the green 
dashed, red dot-dashed and blue dotted curves show results from refs. [TS], [5D] and [51], respectively. 

Since we include the effect of e+e“ decoupling on the neutrino background, described by N^g ~ 3.046 in eq. 
we obtained h [T^g) = 3.9387 in our calculation; here the present photon temperature = 2.7255 ± 0.0006 K 

Our current value of h is thus slightly higher than h (Tly^o) = 3.9138 in ref.[5D] and h = 3.9139 in ref. [15] ■ Note 

that for a given value of the final relic density is directly proportional to 

Because both deconfinement of quarks and gluons and the restoration of the electroweak gauge symmetry are 
associated with smooth cross-overs rather than true phase transitions, the functions g{T) and h{T) are smooth 
everywhere. As noted in the previous Section, we ensured smoothness of these functions by using cubic splines to 
interpolate between different temperature regions. 

There clearly are some differences between the four calculations. These are most visible near the QCD deconfinement 
transition. Moreover, the differences are more visible in gV^, largely due to the derivative term in eq.([^, which 
accentuates the differences between the various treatments. We see that the older calculation used in [TU] 
overestimates gj somewhat for T ~ 0.1 GeV, compared to all three calculations using results from lattice QCD. 
The treatment of ref. [2()j gives a discontinuity at T = Tc, and hence a divergent derivative, yielding a formally infinite 
spike in gj . The continuity has been smoothed out in micrOMEGAs |27j . from which we took the numerical results 
for g and h. We see that this treatment still gives a prominent spike in gl^"^. Apart from this spike, ref.[20] predicts a 
smaller value of in this temperature range than we do. Finally, the prediction for from ref.|21j is quite close 
to our own result, except for some oscillatory behavior just above the QCD transition temperature. 

In order to explore the effect of changes in h and gj on the WIMP relic density, we solve the Boltzmann equation 
1^ numerically. If we choose the initial value Xi of cc to be ^ 10, the final result is independent of the input value 
Y^(xi). We numerically track the behavior of Y^{x) up to a; = 1,000, at which point Y^ has become practically 
constant. The present scaled relic density times squared scaled Hubble constant can then be computed from 


(151, 
42 . 


Pcrit 


f _^^3_ 

SHq/SttGn \ 100 kms“i Mpc“i 


(17) 


where Y^^ = Y-^{x = 1000) and Sq = 2891.2 cm“^; the subscript 0 denotes quantities evaluated at the present time 
[42] . Note that Hq cancels on the right-hand side of eq.(17); this is why the relic density is usually quoted in the 
combination . Numerically, = 2.7889 • GeV). 


1 /2 

The change of the predicted WIMP relic density due to our more refined treatment of the functions h and gj is 
illustrated in fig.^ The upper frame shows results for a temperature independent (cru), while in the lower frame we 
have assumed (cr'iq oc l/x. These behaviors describe the thermally averaged cross section at small velocity away from 
poles (i.e. if the WIMPs cannot annihilate into any particle (p with ~ 2m^) and thresholds (i.e. if the WIMPs 
are significantly heavier than all relevant final-state particles), if the annihilation occurs from a pure 5—wave and 
pure P—wave initial state, respectively. The latter occurs, for example, for a Majorana WIMP annihilating into light 
SM fermions, or for a complex scalar annihilating through s—channel exchange of a gauge boson. Not unexpectedly, 
we observe the largest differences for WIMP masses of a few GeV, which decouple just above the QCD transition 
temperature. The differences amount to up to 9% for the pure S'—wave, and up to 12% for the pure P—wave. 

The results of fig. can be understood in more detail using the approximate analytical solution of the Boltzmann 
equation developed in ref. [15] : 


Vx,o oc -- , XF(x\n (m^gl^^(TF){av){TF)/h(TF)) . (18) 

gY {TF){av){TF) 


1 /2 

Very roughly, xf 20 for WIMP masses and annihilation cross sections of interest. This equation shows that gY 
affects the final result more strongly than h does, which appears only logarithmically. The derivation of eq.(18) 
assumes that gY and h are constant around the WIMP decoupling temperature Tp = m^/xF- This is not a very 
good approximation near the QCD deconfinement transition, where these functions change rapidly, as we saw in fig.[l] 
However, we can see directly from the Boltzmann equation that the most relevant temperature range is that around 
the decoupling temperature. At higher temperatures, Y^ is in any case close to its equilibrium value, which does not 
depend on gY ■ At temperatures well below the decoupling temperature, i.e. for x ^ xf, the right-hand side of the 
Boltzmann equation ([^ is suppressed by the explicit x~'^ factor. If {av) oc l/cc, as in the lower frame of fig. the 
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FIG. 2: The relative difference between the predicted relic density of a Majorana WIMP between our calculation and a 
calculation using the same older results for the functions h and shown in Fig. as function of the WIMP mass. The 
upper frame is for a constant (cru), chosen such that our prediction for = 0.1193, while the lower frame is for a pure 

P—wave annihilation, with {av) = 1.2 • 10“^'* cm®s“^ • Tjra^. These results are almost independent of the numerical size of 
the annihilation cross section. 


1/2 

suppression at a; > a;i? is even stronger. Sharp features in gj therefore give sharper features, with larger amplitudes, 
for pure P—wave annihilation than for 5'—wave annihilation. 


We noticed earlier that the older treatment of ref. m overestimates for some range of temperatures above Tc- 
Eq.(18l indicates that this should lead to a smaller predicted relic density, which is confirmed by fig. Since gV^ is 
over-estimated for an extended range of temperatures, the effect on the relic density is about the same for S— and 
P—wave annihilation, amounting to about 5% near the peak of the ratio shown in fig. The second, much lower 
peak near = 1 TeV is probably due to lack of knowledge of the top mass at the time when ref. [H] was written. 

The spike in gj predicted by the micrOMEGAs treatment of the results of ref. [20] gives prominent spikes in the ratios 
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< 20 GeV. Finally, we note that treating the charm quark as a free 


shown in figj^ These spikes are numerical artefacts that result from the smoothing procedure used in micrOMEGAs. As 
argued in ref. [20] , a true S function spike in should not affect the numerical solution of the Boltzmann equation, 
which necessarily entails some discretization. The probability that the program then has to evaluate the right-hand 
side of the Boltzmann equation at the precise value of x where the S function diverges is zero. We nevertheless 
show results including this spike since it results from the “standard treatment” encoded in micrOMEGAs. Outside the 
mass range affected by this spike, the results of ref. [20j predict a slightly too large relic density, consistent with our 
observation that it predicts smaller values of gj and h than our treatment does. Note that for fixed gj , reducing 
h will (slightly) increase xp, leading to an increase of the predicted relic density. A decrease of h therefore goes into 
the same direction as a decrease of gj . However, the fact that the relative difference between our calculation and 
the prediction based on ref. [20] is almost the same in both frames of fig. at large WIMP masses shows that the main 
effect still comes from the change of gl^"^■ 

We saw in fig. that the prediction for g^^ from ref. m lies below our prediction, except for a very narrow range 
of temperatures around Tc- As a result, for pure S'—wave annihilation the prediction for the relic density based on 
the treatment of ref. [21] lies above our prediction for all WIMP masses larger than 2 GeV. We argued above that the 
relevant range of temperatures is (even) smaller for pure P—wave annihilation. This explains why the blue curve in 
the lower frame of fig. goes slightly above 1 for ~ 25 GeV. Note also that the predictions using our treatment 
agrees with the prediction using ref.[2T] to better than 1% for all WIMP masses, except in the range between 3 and 
15 GeV where the difference reaches 9 (12)% for pure S — (P—)wave annihilation. 

In order to put these results into perspective, it should be noted that the lattice QCD predictions for the energy 
and entropy densities listed in Table 1 of [23] . on which our treatment is based, still have significant uncertainties, 
which decrease from about 14% at T = 130 MeV to about 3% at P = 400 MeV. The corresponding uncertainty in 
the relic density is up to 2.5% for 2 GeV < m 

particle would increase the predicted relic density by about 2.2% for ~ 30 GeV. 

The results shown in fig. [^ are almost independent of the assumed value of the WIMP annihilation cross section as 
long as the relic density comes out at least roug hly correctly, although they evidently are sensitive to the functional 
dependence of (av) on the temperature. Eq.(18l shows that thermodynamic effects enter primarily through gj {xp), 
and xp depends on the annihilation cross section only logarithmically. 

On the other hand, the exact value of the annihilation cross section that reproduces the correct relic density, now 
(within standard ACDM cosmology) constrained to be = 0.1193±0.0014 [331: does depend on gV^{T ~ Tp) and, 
to a lesser extent, on h{T ~ Tp). Precise knowledge of this cross section is important to constrain the free parameters 
of models of thermal WIMPs. Moreover, as we will see in more detail below, indirect DM searches now begin to 
probe annihilation cross sections close to the required value of {crv), if the latter is (approximately) independent of 
the temperature. 

In fig. we show the required value of (crv), assumed to be independent of the temperature, for a Majorana fermion, 
obtained from our refined calculation of gl^^ and h. This updates the results of ref.[33]) which assumed = 0.11 
and used [H] to compute gl^"^ and h. 

We see that for 10 TeV > > 10 GeV the required value of {av) is in fact closer to 2 • 10“^® cm^s“^ than to the 

“canonical”, often cited value of 3 • 10“^® cm®s“^. The near constancy of the required value over such a large range 
of WIMP masses results from an “accidental” cancellation of two effects. This can again be understood from the 
approximate analytical solution (18) of the Boltzmann equation. On the one hand, increasing increases xp, which 
increases the relic density. Since xp depends only logarithmically on m^, the freeze-out temperature Tp = m^/xp 

still increases as is increased. As shown in fig. 1 this increases gV^{Tp), which in turn reduces the relic density. 

^ 1/2 

For > 10 TeV, all SM particles are essentially fully relativistic at Tp, i.e. gT becomes independent of T, reaching 
its asymptotic value of 106.75. For these very large WIMP masses the required value of {av) would thus increase 
logarithmically with in order to cancel the effect of the increase of a;i?. However, since by dimensional analysis 
and unitarity arguments [33] {av) oc l/m^, it is very difficult to find scenarios with sufficiently large WIMP mass for 
> 10 TeV. 

On the other hand, for WIMP masses below 10 GeV the rapid decrease of gV^{Tp) with decreasing Tp shown in 


fig. [B requires a rather rapid increase of {av), to a peak value of about 4.5 • 10 cm®s Finally, for < 0.35 

GeV, gj {Tp) becomes approximately constant again, with electrons, positrons, neutrinos, and photons contributing 
1 /2 

so that gT ~ 3.29. Since xp keeps decreasing with decreasing m^, keeping the relic density constant requires that 
{av) also decreases logarithmically with decreasing WIMP mass for these very light WIMPs. 
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FIG. 3: The value of {(Jv), assumed to be completely independent of temperature, required to obtain a thermal relic density 
= 0.1193 within standard cosmology, as a fnnction of WIMP mass. 


4. EXPERIMENTAL CONSTRAINTS ON (av) 

In this section, we will compare experimental constraints from indirect WIMP searches and from analyses of the 
cosmic microwave background (CMB) with our prediction for (av) shown in fig. Since the CMB decoupled much 
later than WIMPs did, and hence also at a much lower temperature 0.3 eV rather than ~ m^/20), while WIMPs 
in galaxies now have an average kinetic energy of ^ such a comparison is meaningful only if (av) is largely 

independent of the temperature. If {av) oc T, as in pure P—wave annihilation, or for even stronger T—dependence, 
the bounds on {av) from the CMB and from indirect WIMP searches are still several orders of magnitude above the 
value required to obtain the correct relic density. 

Currently the strongest and most robust upper bounds on {av) from indirect WIMP searches come from searches 
for hard 7 rays by the FermiLAT collaboration. Photons travel in straight lines through our galaxy, whereas charged 
particles get deflected by the galactic magnetic field. This not only introduces a sizable uncertainty, since our 
knowledge of this magnetic held is far from perfect; it also isotropizes the arrival directions of the WIMP annihilation 
products, making it impossible to focus on regions of space where the WIMP signal should be particularly strong. 
Another advantage of 7 rays is that they are present in nearly all possible hnal states: hard photons can be emitted 
directly off hnal or intermediate state particles, can originate from the decay of neutral pions and other hadrons that 
result from the hadronization of qq hnal states or the decay of r leptons, and can be produced from energetic electrons 
or positrons through “inverse Compton” upscattering of ambient photons. 

The strongest WIMP signal is expected from near the center of our own galaxy. Unfortunately this region also 
hosts several backgrounds, both in form of point sources and in form of extended emission. It has been claimed that 
there is evidence for an additional component in the GeV 7 hux from near the galactic center which can be explained 
through WIMP annihilation |32j : but other interpretations of this additional component exist [15] . We also note that 
the FermiLAT collaboration itself has not published any analysis of their data on the galactic center. 

In this paper we therefore focus on FermiLAT observations of nearby dwarf galaxies [T3HT6] . In contrast to big 
galaxies like our own, the mass density of dwarf galaxies should be dominated by dark matter even in the central 
region, yielding a much better signal-to-background ratio for indirect WIMP signals. No such signal has been seen. 
Our analysis is based on the very recent 6 -year “Pass 8 ” analysis m- 

The results are shown in fig. We see that the upper bound on {av) is strongest if WIMPs predominantly 
annihilate into uu final states, but the bound for WIMP annihilation into bb is only slightly weaker. For the r+r” 
final state the upper bound on the cross section is similar for WIMP masses below 40 GeV, but is somewhat weaker 
for heavier WIMPs; hadronic final states have higher multiplicity, and hence higher 7 flux per WIMP annihilation, for 
larger WIMP masses, whereas for the T+r“ final state the photon multiplicity is essentially independent of the WIMP 
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mass. These constraints exclude WIMPs with mass < 70 to 100 GeV annihilating into hadrons or r leptons with 
temperature independent (erf). 



FIG. 4: The result of fig.is compared with several observational upper bounds on (cru), which is assumed to be independent 
of temperature. The cyan, green and magenta curves follow from the FermiLAT upper bound m on the 7 flux from dwarf 
galaxies, for different dominant WIMP annihilation channel (uw, bb or r'*’T“), whereas the red curve results from an upper 
bound on spectral distortions of the CMB, assuming WIMP annihilation into e^e“ pairs. 

As mentioned above, WIMP annihilation into e“'"e“ can produce hard photons through upscattering ambient, e.g. 
visible, photons. The corresponding upper bound on (av) (not shown) is worse than that for WIMP annihilation into 
T+T“ by a factor of about two to three m, excluding WIMPs with mass < 15 GeV annihilating into e“''e pairs 
for the value of (av) shown in fig. 

WIMP annihilation can also affect the CMB. The strongest limits on (av) originate from the fact that WIMP 
annihilation heats up the plasma in the “recombination” epoch when neutral atoms first formed m, thereby delaying 
the decoupling of the CMB photons and distorting the pattern of CMB anisotropies. In fig. we show the bound on 
the WIMP annihilation cross section into e“'"e“ pairs that results from an analysis m of data from the WMAP and 
ACT collaborations. It excludes a thermal WIMP with < 8 GeV. 

PLANCK data will lead to considerable stronger constraints [131I1S]. Unfortunately these papers only cite upper 
bounds on the product of the WIMP annihilation cross section and an efficiency factor fes with which the energy of 
the WIMP annihilation products is absorbed in the thermal plasma. Using results from ref. [H] , we estimate that the 
latest PLANCK data exclude WIMPs with 40 GeV annihilating into e“''e“ pairs with temperature independent 
cross section; see also the recent analysis [50] . 

Since the efficiency factor should be similar for the other final states considered in fig. the CMB constraint should 
also be similar for all channels. The current CMB constraint is thus weaker than the bounds derived from the most 
recent FermiLAT data if WIMPs mostly annihilate into qq or T+r“ final states, but is stronger for WIMPs annihilating 
predominantly into e"'"e“. However, one should keep in mind that the CMB constraint is less direct. It is conceivable 
that additional non-standard ingredients to the CMB fit - e.g., the presence of sterile neutrinos, a significant running 
of the spectral index of inflation, and/or a large contribution from tensor modes - can (partly) compensate the 
distortions caused by early WIMP annihilation, thereby weakening the constraint on {av). On the other hand, the 
constraint derived from the observation of dwarf galaxies depends on the assumed dark matter distribution |13| . In 
any case, it is encouraging that recent astrophysical and cosmological observations begin to probe relatively light 
thermal WIMPs with temperature independent annihilation cross section. 
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5. SUMMARY AND CONCLUSIONS 

Using recent lattice QCD results with dynamical quarks for the equation of state, we have computed the energy 
and entropy densities of the SM with emphasis on temperatures around the deconfinement transition at Tc = 154 
MeV. These results are described by the functions g{T) and h{T). Of particular relevance for the calulcation of the 
relic density of thermal WIMPs is the quantity defined in eq.(|^, which also depends on the derivative of h with 
respect to the temperature. Our results for these functions can readily be embedded in the public codes computing 
the relic density [26U28] . 

Our predictions for the WIMP relic density differ from earlier treatments that relied on phenomenological models 
or pure glue lattice QCD calculations. These differences are most pronounced for WIMP masses between 3 and 15 
GeV; they can reach about 9% for a constant (av), and up to 12% if (av) oc T as in pure P—wave annihilation away 
from poles and thresholds. These differences are partly due to a spurious divergence in h{T) leading to a divergence in 
g*' , which results in a sharp spike in the standard implementation. This illustrates the importance of ensuring that 
h{T) is continuous everywhere, not just during the deconfinement transition but also during electroweak symmetry 
breaking at T ~ 100 GeV, which we essentially ignored in our treatment. 

It should be noted that the uncertainties in the recent lattice calculation we used still translates into an uncertainty 
of the predicted relic density of up to 2.5%, considerably larger than the observational uncertainty (at least in the 
framework of the standard ACDM cosmology). On the other hand, for > 20 GeV our result for the relic density 
agrees to better than 1 % with the best previous calculation pT] . 

We used our improved treatment of the thermodynamics of the very early universe to update the calculation of 
the WIMP annihilation cross section required to reproduce the observed dark matter relic density, assuming the 
thermal average {crv) to be independent of temperature. We found that this cross section is indeed nearly constant 
for 10 GeV < < 10 TeV, thanks to a fortuitous cancellation between two competing effects; however, as also 

pointed out in ref. [24], the required value is closer to 2 • 10“^® cm^s“^ than to the often-cited “canonical” value 
of 3 • 10“^® cm^s“^. On the other hand, for ^ 3 GeV the required value of (av) exceeds 4 • 10“^® cm^s“^. 
Producing a sufficiently large annihilation cross section for such light WIMPs typically requires the introduction of 
new light “mediators” between the WIMPs and SM particles. If the mediator mass is ^ m^/10, the contribution of 
the mediator particles to the entropy and energy densities should be included, slightly modifying our results for the 
required annihilation cross section shown in fig. [^ 

We also compared the required value of {av) with upper bounds on this quantity that come from searches for 
energetic 7 rays produced in WIMP annihilation as well as from CMB constraints. The most recent PLANCK data 
most likely give the strongest upper bounds, although a dedicated analysis for specific final states has not yet been 
performed. Moreover, the CMB bound assumes standard ACDM cosmology, and can thus more easily be evaded 
than the bounds from 7 ray searches. All of these bounds can be evaded if the WIMPs predominantly annihilate into 
neutrinos m- 

It should be noted that both our calculation of the required WIMP annihilation cross section and our analysis of 
observational upper bounds on this quantity assumed that the thermal average (av) is independent of temperature, or, 
equivalently, that the annihilation cross section is independent of the invariant center-of-mass energy. Theoretically 
this assumption is not well motivated. For non-relativistic WIMPs the annihilation cross section can usually (but not 
always [52] '! be expanded in terms of the relative velocity v, av = a + bv^ + • • •. Even if the constant term a is not 
suppressed, i.e. if annihilation from an S'—wave is allowed, one would generically expect b to be of the same order 
as a. This would reduce the required value of (av) in today’s universe by about 10%. If the constant term in av 
is suppressed, the upper bounds on (av) that follow from observations in today’s universe do not constrain thermal 
WIMP models significantly. 

However, even in this case it is important to calculate the relic density as accurately as possible. This leads to a 
constraint on the free parameters of the underlying WIMP model, which can hopefully one day be compared to direct 
measurements at colliders. Only then will we be able to say with some confidence that this WIMP model is indeed 
correct. This paper makes a small contribution to this ambitious long-term program. 
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